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We model a radiating, moving black hole in terms of a worldtube-nullcone boundary value prob- 
lem. We evolve this data in the region interior to the worldtube but exterior to a trapped surface by 
means of a characteristic evolution based upon a family of ingoing null hypersurfaces. Data on the 
worldtube is induced from a Schwarzschild spacetime but the worldtube is allowed to move relative 
to the static Schwarzschild trajectories. When the worldtube is stationary (static or rotating in 
place), a distorted black hole inside it evolves to equilibrium with the Schwarzschild boundary. A 
boost of the worldtube with respect to the Schwarzschild black hole does not affect these results. 
The code also stably tracks an unlimited number of orbits when the worldtube wobbles periodically. 
The work establishes that characteristic evolution can evolve a spacetime with a distorted black hole 
moving on a 3-dimensional grid with the controlled accuracy and long term stability necessary to 

ON , 

investigate new facets of black hole physics. 
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r—i ■ 

q ■ I. INTRODUCTION 

o : 

The calculation of gravitational waveforms from the inspiral and merger of binary black holes is a great current 
challenge to computational relativity with important relevance to LIGO. A three dimensional computer code is now 
being constructed by the Binary Black Hole Grand Challenge Alliance to solve this problem by Cauchy evolution of 
, initial data for two black holes pj . The most difficult problems with the development of this code are inaccuracy and 
^ ■ instability at the boundaries of the Cauchy grid. Matching to an exterior characteristic evolution is one method being 
pursued to handle the outer Cauchy boundary and to extract the waveform at null infinity. This has been shown to 
be a highly accurate and efficient approach in the treatment of 3-dimensional nonlinear scalar waves For the 

purpose of extending this approach to full general relativity, a 3-dimensional characteristic code for the gravitational 
field has been developed and fully calibrated to perform with second order accuracy and robust stability in highly 
j-^. nonlinear regimes A module for matching Cauchy and characteristic gravitational evolution codes across a 

Q\ , worldtube interface has also been written |Q and is now in the testing stage. In this paper, we present results which 
show that Cauchy-characteristic matching can also solve the very difficult inner boundary condition necessary for 
q-( Cauchy evolution of black holes. Notably, we have used a characteristic code to achieve the first successful treatment 
I ■ of a distorted black hole moving on a 3-dimensional grid with apparently unlimited long term stability. 

The conventional strategy for avoiding the topological and strong field difficulties in the Cauchy evolution of black 
holes has been to excise the region interior to an apparent horizon, as initially suggested by W. Unruh |H|. In a recent 
work j^] , we proposed an alternative version of this strategy in which the black hole region is evolved by a characteristic 
evolution based upon ingoing null cones. These null cones are truncated at an inner boundary consisting of a trapped 
or marginally trapped surface, and matched at their outer boundary to the inner boundary of a Cauchy evolution. 
In turn, the outer boundary of the Cauchy evolution is matched to an exterior characteristic evolution extending to 
(compactified) infinity. The potential advantages over a purely Cauchy approach to the inner boundary were discussed 
in that work and the global strategy was successfully implemented for spherically symmetric self-gravitating scalar 
waves evolving in a black hole spacetime. 

We present here the successful implementation of a characteristic treatment of an inner black hole boundary for 
fully 3-dimensional simulations containing gravitational radiation. We show that the ingoing characteristic approach 
is able to locate the black hole and to track it stably as it moves on the numerical grid. For a report on recent 
progress in tackling the same problem by Cauchy evolution see Ref. ]To| ]. The boundary data for the characteristic 
initial value problem is posed on a worldtube and on an ingoing null cone emanating from the initial slice of the 
worldtube. A main goal of this study is to develop new methods which will allow a combination of characteristic and 
Cauchy evolution to tackle the computational problem of the inspiral of two black holes. In this new approach, the two 
evolutions are matched across a worldtube, with the Cauchy domain supplying the boundary values for characteristic 
evolution and vice versa. In treating this problem, there are major computational advantages in posing the Cauchy 
evolution in a frame which is co-rotating with the orbiting black holes. Indeed, such a description may be necessary in 
order to keep the numerical grid from being intrinsically twisted and strangled. In this co-orbiting formulation of the 
binary black hole problem, the Cauchy evolution requires inner boundary conditions in two regions approximating the 
two disjoint apparent horizons and also an outer boundary condition on a worldtube. Far enough outside the outer 
worldtube the coordinate rotation would go superluminal. Previous work has shown that an outgoing characteristic 
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code can routinely handle such superluminal gauge flows in the exterior. Our present results indicate that an ingoing 
characteristic code can just as effectively handle the inner boundaries of multiple black holes. 

Consistent worldtube data must satisfy conservation conditions j||5) which correspond to a hyperbolic version of 
the standard constraint equations for Cauchy data. In the present work, since we are not matching to a Cauchy 
evolution, we generate this data from an analytic solution. Specifically, we set Schwarzschild data at the worldtube. 
In ingoing Kerr-Schild coordinates x a , the Schwarzschild metric takes the form || 

,? ,72 >~2 ,-2 ,~? 2m. - xdx + ydy + zdz , 2 . 

ds 2 = -dt 2 + dx 2 + dy 2 + dz 2 + —(dt+ ^ ) , (1.1) 

r r 

where f 2 = x 2 + y 2 + z 2 and fc M = —d^(i + f) is the ingoing degenerate null vector. 

In order to determine data for a moving black hole, we introduce an x a coordinate system which is either rotating, 
boosted or wobbling with respect to the x a coordinates. In the x a frame, the null coordinates are centered at r = 
and the worldtube is located at r = R, where r 2 = x 2 + y 2 + z 2 . As a result, the worldtube is stationary with respect 
to the null grid but the black hole moves. 

The initial value problem also requires data on the initial null hypersurface. In our formulation, this null data 
is free of constraints, other than continuity conditions at the worldtube. Thus we can introduce an arbitrary pulse 
of radiation in the data to describe a distorted Schwarzschild black hole. We can also pose initial null data by 
setting to zero the component of the Weyl tensor intrinsic to the null hypersurface. On a spherically symmetric null 
hypersurface, the Weyl data for the Schwarzschild space-time is exactly zero. On a null hypersurface offset from the 
center of spherical symmetry, this Weyl data for a Schwarzschild spacetime is not zero (and it is not possible to express 
it in simple analytical form). In this case the choice of vanishing Weyl data introduces an amount of gravitational 
radiation which depends on the size of the offset. The details of setting up the worldtube- nullcone data are presented 
in Section O. 

The worldtube boundary values, as prescribed in this paper, satisfy the conservation conditions when the spacetime 
is exactly Schwarzschild, e.g. a Schwarzschild black hole in either rotating or boosted coordinates. But in either the 
distorted or wobbling case, when gravitational radiation is contained in the initial null data, the extraction module 
overdetermines the metric and its normal derivatives at the worldtube in terms of their Schwarzschild values. As a 
result, the reflection of the wave off the worldtube can lead to a violation of the Bianchi identities. The mismatch 
between radiation impinging on the Schwarzschild worldtube introduces an unphysical sheet source of gravitational 
radiation on the worldtube which is not necessarily consistent with energy and angular momentum conservation. 
Although this obscures the physical interpretation of the results for those cases, it is remarkable that the stability 
of the evolution is not affected and that the system behaves in accord with the principles of black hole dynamics, 
as described in Sec. [v|. In a more physical implementation, the conservation conditions would be enforced either (i) 
directly, by using them to properly evolve the boundary conditions up the worldtube, or (ii) by matching the interior 
evolution across the worldtube to an exterior Cauchy evolution. 

Our chief goal here is to demonstrate that the region of spacetime interior to the world tube can be evolved stably 
and accurately by means of a characteristic evolution algorithm. The interior region is truncated in the vicinity of 
an apparent horizon. The long term stability we observe indicates a surprising robustness of the worldtubc-nullconc 
boundary value problem. 

In the case of a Schwarzschild spacetime described in a frame rotating about the center of spherical symmetry, the 
location of the apparent horizon is known analytically, as well as the transformation to null coordinates and the null 
metric. Thus this case provides an important test bed to calibrate numerical accuracy. Long term stability and second 
order convergence to the analytic values have been confirmed. In this purely rotating case in which the worldtube is a 
stationary boundary, when we superimpose a pulse of radiation on the initial Schwarzschild null data we find that the 
surface area of the resulting distorted black hole grows in time but eventually reaches an equilibrium value consistent 
with the Schwarzschild boundary conditions on the worldtube. In the offset case, the Schwarzschild boundary moves 
periodically but the marginally trapped surface associated with the black hole again reaches equilibrium with it, 
confirming that the motion of the boundary is "pure gauge" . 

When the null cones are not spherically symmetric a computational approach is necessary to find the "trapping 
boundary" , where a marginally trapped surface is located on each ingoing null hypersurface. The analogous problem 
in Cauchy evolution is the excision of a black hole interior by locating the apparent horizon. There is an extensive 
literature on marginally trapped surface (MTS) finders on spacelike hypersurfaces |TT|-po||. In Section III, we present 



two methods for use on null hypersurfaces. Computational design and performance is discussed in Sec. IV. 
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II. DETERMINATION OF THE DATA 



The worldtube data fo r th e null evolution of the interior region is obtained by recasting the Cartesian form of the 
Kerr-Schild metric Eq. (Ol) for a Schwarzschild black hole into nonstationary ingoing null coordinates. However, 
the Schwarzschild metric can only be expressed in an analytically manageable null coordinate form when the null 
cones are centered to be spherically symmetric. Thus, in the offset case, numerical techniques must be used to carry 
out the transformation from Cartesian to null coordinates in order to provide the worldtube data. Fortunately, the 
transformation need only be implemented in the neighborhood of the world tube, This allows the solution for the 
ingoing null geodesies to be carried out by means of a Taylor expansion in affine parameter, up to the required order 
of accuracy. In this way, the transformation can be formulated as a hybrid analytical-numerical scheme, with the null 
properties built in analytically. 

Such a scheme for extracting null worldtube data from Cauchy data in a general spacetime has been implemented as 
an extraction module [Q , which is part of the computational procedure to match Cauchy and characteristic evolutions 
across an interface. We use this module here to obtain the required worldtube data for null evolution by extraction 
from the "3+1" form of the Schwarzschild metric which has been recast into nonstationary coordinates by a nontrivial 
choice of lapse and shift . 

The initial value problem is completed by giving data on an initial null hypersurface. In ingoing null coordinates, 
with v labeling the null hypersurfaces, with x A — (8, <f>) labeling the angles along the null rays and with r labeling 
the surface area measured along the null rays, the metric takes the Bondi-Sachs form |)| 

ds 2 = (e 2/3 ^ + r 2 h AB U A U B ^j dv 2 + 2e 2f3 dvdr - 2r 2 h AB U B dvdx A + r 2 h AB dx A dx B , (2.1) 

where det(h,AB) = det{q AB ) = <7, with q AB a unit sphere metric. The inverses of these 2-dimcnsional metrics 
are denoted by h AB and q AB . We express q A B in terms of a complex dyad q A (satisfying q A q A = 0, q A q A = 2, 
q A = q AB qB, with q AB q B c = 5 A and q A B = \ {qAQB + QaIb))- h A s can then be represented by its dyad component 
J = h A Bq A q B /2, with the spherically symmetric case characterized by J = 0. The full nonlinear h A s is uniquely 
determined by J, since the determinant condition implies that the remaining dyad component K = h A sq A q B /2 
satisfies 1 = K 2 — J J. We also introduce the spin- weighted field U = U A q A , as well as the (complex differential) eth 
operators 3 and 8. Refer to j|,^2| for further details. 

The complete null data can be specified freely in terms of either the metric quantity J or the Weyl tensor component 
C a f3 1 sn a n' y m^m s (corresponding to '5 4 in Newman-Penrose terminology |23[| ), where the null vector n a and the 
complex spacelike vector m a span the tangent space to the null hypersurface. On a spherically symmetric null 
hypersurface, the Weyl data for the Schwarzschild space-time is exactly zero but since there are no constraints on the 
null data, we can freely add a radiation pulse of any desired shape. 



In the case of rotating coordinates, as described in Sec. |V A| , or boosted coordinates, as described in Sec. |VB] , 
the initial null cone is spherically symmetric and Schwarzschild null data can be determined analytically. However, 
it is not possible to express Schwarzschild null data in simple analytical form on a null hypersurface which is not 
spherically symmetric. Instead, we choose initial data for an offset wobbling black hole by setting the Weyl data to 
zero on the initial (nonsymmetric) null hypersurface. The relevant components of the Riemann tensor are 

tirABr — -^9 AB .rr ~ ^9 9vr,r9AB,r ~ ~^9 9AC,r9BD,r 

= -^h AB ,rr + rh AB ^ r - r(3 <r (rh AB<r + 2h AB ) - —h CD h AC , r h BD ^. (2.2) 



Here g AB R rABr = by virtue of the hypersurface equation for (3 0, 



Ar = -^rh AB h ABtr = ^(J, r Jr ~ K f r ) . (2.3) 

The requirement of vanishing Weyl data is equivalent to q A q B R rABr — 0, which gives 

r 2 

{r J, r ),r ~ 2/3, r (r 2 J), r - — j(j r J, r - K 2 r ) = 0. (2.4) 

Combining Eq's. ( |2.3| ) and (p~i|), we then obtain 

r 2 (r 2 J r ), r - 2/3 r (r 4 J), r = 0. (2.5) 
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With /3 r determined from Eq. ( |2.3| ), Eq. (2.5) is a second order radial differential equation for the initial data J, 
which may be solved in terms of boundary values for J and J r on the worldtube. 

The outer worldtube is located at (x 2 + y 2 + z 2 ) 1/>2 = R, in terms of coordinates x a = (t,x,y,z) moving with 
respect to stat ic Kerr-Schild coordinates x a . The boundary data on the worldtube is obtained by first transforming 



the metric (1.1) in "3+1" form to the x a frame and then applying the extraction module, which determines numerically 
the transformation from a Cartesian to a null coordinate system in the neighborhood of the worldtube. This module 
supplies the boundary values of the null metric quantities J, (3, U and V on the worldtube. 

As a check on the extraction module, we examine the rotating case where the null Schwarzschild metric can be found 
analytically. We relate the static coordinate frame x a to the rotating one, x a by t — i , x — xcosoj(t + r) — y smui(t + r), 
y = x sinu)(t + r) + ycosu(t + r) and z = z. In this transformation, the angles 9, 4> (associated in the standard way 
with the Cartesian coordinates x l ) remain constant along the generators of the null cones emanating from the world 
tube. Therefore, the null metric can be easily obtained by the simple transformation v = t + r, r = \J (x 2 + y 2 + z 2 ), 
6 = cos _1 (z/r) and <f> — taxi~ 1 (y/x), leading to 

ds 2 = ( - 1 + + r 2 uj 2 sin 2 9 \dv 2 + 2dvdr - 2t 2 uj 2 sin 2 9dudc/> + r 2 (d9 2 + sin 2 9d(f) 2 ). (2.6) 



This identifies the spin weighted versions of the variables appearing in the null metric (2.1) as 



J = 0,/3 = (2.7) 
U = i lu sin 9 e l4 (2.8) 
V = -r + 2m. (2.9) 

We used this transformation to check the accuracy of the worldtube data extracted from the Schwarzschild metric 
in "3+1" Cauchy form with mass m = .25. The code was run using the extraction radius R = 3. The Cauchy grid 
was a Cartesian cube with range x % £ [—4,4] and the range of the characteristic grid was r £ [0,4]. We confirmed 
convergence of the numerical error to zero at a second order rate with grid size for several values of u> (0.8, 0.1, 0.001). 

III. FINDING THE TRAPPING BOUNDARY 

The excision of a region inside the black hole is necessary for numerical evolution, for otherwise contact with the 
singular region of spacetime would result. The boundary of the excised region must either be trapped or marginally 
trapped in order to insure that it does not contain points that can causally influence the radiated waveform. For 
a slice S of an ingoing null hypersurface J\f v , described in null coordinates by r = R(v,x A ) 7 the divergence of the 
outgoing null normals is || 

= -V- -±-[^{e 2 ^h AB R^ B - r 2 U A )U - r(r~ 1 e 2 Pfi AB ) jr R t AR,B + r 2 U A R, A . (3.1) 
This is to be evaluated on S after all partial derivatives are taken; e.g. r a — 0. The slice will be marginally trapped 

if e ; = o. 

Finding a marginally trapped slice on a convergin g ing oing null hypersurface is a 2 dimensional elliptic problem, 



which entails setting the right hand side of equation (3.1) to and solving for R(v,x ), at a fixed advanced time v 



It is easier to find trapped surfaces. In fact, the largest r — const slice of Af v that satisfies the algebraic inequality 
Q < 0, where 

Q = -V + ^={y/qU A \A. (3.2) 

i s ci ther trapped or marginally trapped (9). We call this slice the "Q-boundary" . A comparison of Eq's. ( |3~l| ) and 
( |3.2|) shows that the Q-boundary is marginally trapped when the slice Q — is an r = const slice. However, the 
gauge freedom in the choice of a surface area coordinate (r is a scalar density) allows any slice to be regauged as an 
r = const slice. So there is a gauge in which the Q-boundary and the trapping boundary coincide. But finding this 
gauge is tantamount to solving the elliptic problem for a marginally trapped slice. 

This presents us with two possible strategies for locating the inner boundary, both of which ensure that the excised 
portion of spacetime cannot causally effect the exterior spacetime: (I) Use the trapping boundary or (II) use the 
Q-boundary. Strategy (I) makes the most efficient use of the spacetime points but a 2D elliptic equation must be 
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solved. Strategy (II) involves only some simple algebra so it is very efficient computationally. We will pursue both 
strategies here and compare their merits. 

In implementing strategy (I), we have tried two methods for finding the trapping boundary. One is a standard 
approach to the solution of the elliptic equation 0; = by solving the parabolic equation 

d x R = -Qi (3.3) 

in terms of a relaxation parameter A. At large A, the solution relaxes to the location of the trapping boundary 
provided that the procedure is stable. 

Another method of finding the trapping boundary is by a minimization procedure. In the case of a Cauchy 
hypersurface, this approach was introduced in [p"8"[and further developed in [ p^[ using spectral methods. Here we use 
a finite difference version of the minimization approach and combine it with an approach based upon gradient flows 
proposed in |l3( ] . We combine these approaches by characterizing the trapping boundary on a converging ingoing null 
cone Af v as a marginally trapped slice S m which minimizes the functional 

T(S) = l<f e 2 n efds (3.4) 



where O ra is the divergence of the ingoing null vector n a tangent to M v . We normalize n a and l a (both normal to S) 
by n a l a — — 1. With this setup, T is invariant under changes in the extensions of the the null normals (local boosts) 
that retain n a l a = — 1. For a standard sphere of radius R in Minkowski space, T — 8ir/R 2 so that the minimum on 
a Minkowski light cone occurs at infinity. This has the advantage of biasing the search away from the caustic tip of 
the light cone when looking for nontrivial minima. 

We perform a variation of the form 5x a = Fn a 5X, which deforms the slice along the generators of J\f v . In order to 
find a flow F which leads toward the minimum consider the variation 

8T(S) = j> T a Sx a dS, (3.5) 

which serves to define T a . We choose F = —T a n a . Then 

ST(S) =-j> {T a n a ) 2 SXdS, (3.6) 

so that the variation will lead in the direction of the minimum at T = 0. 

The main problem then reduces to calculating T a n a . In ingoing null coordinates on Af v , we describe S by r = R(x ) 
and its variation by r = R(x A ) + SR(x A ), with Sv — Sx A = 0. Choosing the extension n a — — g vr v ia , we have 
n = — 2/r so that, for any slice of J\f v , 8Q n = 2Sr/r 2 . 

The variation of terms f(r,x A ) not explicitly dependent on R is calculated using 5f(r,x A ) — f^SR- Also, since 
dS = r 2 df2, in terms of the solid angle element in the x coordina tes, we have 5dS — 2dSSR/r. As a result, the 



contributions from 5Q n and from SdS cancel in the variation of Eq. (3.5) so that 



ST(S) = j> BiSBidn, (3.7) 



From equation ( 3T ) 



Op-2/3 1 

SQt = Qi.rSR + -^{ -Wq{e 2fi h AB 8R,BlA - 2r{r- 1 e 2(3 h AB ). r R. A SR, B + r 2 U A 5R A ] (3.8) 

r v9 

For any vector field V (r, x B ) on S, we have 

D A V A dn = 0, (3.9) 



is 

where DaV a = (^/qV A ),A/\/q + V a R,a- This allows us to eliminate terms in ST containing SR } a by integrating over 
S. We obtain 

8T(S) = j> {@i@i, r + ^jSRdn (3.10) 
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where 



* = [jsCejO.r " ^{Qie- W ),re^h AB l r R, A R,B - [2@ie-^U A ], r R, A - {2@ l e-^U A ], A 

l^(®ie- 2P ),Ae 20 h^ ir R B + [(^Q lh AB ), r R A - (% e -^), r ^-h AB R <A - ^(O^^h^U- (3.11) 
Thus, in order to find the trapping boundary, we follow the variational path determined by SR = —FSX with 

F = (QiQi, r + ^)/R 2 . (3.12) 

As a check on the stability of this minimization scheme, suppose iS m is a trapping boundary located at r = R m (x A ). 
Then, on S m , 0/ = and F = ty/Ri^. But the operator (d A + R m , A d r ) also annihilates 9; on <S m . As a result, direct 



substitution to eliminate the d r derivatives in equation (3.11) gives \& = on S m . Thus 5R also vanishes and S m is 
an invariant slice with respect to the variational scheme. 
We have: 

(1) A marginally trapped surface is a zero of the positive functional T; 

(2) The effect of the flow F on T is everywhere negative or zero; 

(3) A marginally trapped surface is stationary under the flow F. 

Thus a marginally trapped surface must be stable with respect to the flow F except in the degenerate case, 
corresponding to neutral stability, where a continuum of such surfaces exist in Af v . Although this degenerate case 
is possible it would be improbable to encounter in an evolution based upon a reasonably behaved foliation. For an 
interesting discussion of the wild behavior possible in general for marginally trapped surfaces see |H| . 

In order to implement either of the above two finders as computational algorithms, we represent the geometric 
quantities involved as spin-weighted fields in stereographic coordinates. The spin-weighted expressions necessary to 
determine Oj and F are given in Appendix [a|. 



IV. COMPUTATIONAL DESIGN AND PERFORMANCE 



All numerical algorithms have been based upon explicit finite difference methods. The metric functions are dis- 
cretized and placed on a finite 3 dimensional grid, with N r radial points and 7V| angular points, whose outer boundary 
is a spherical worldtube. The spherical coordinates are patched by two overlapping stereographic grids and angular 
derivatives of tensors are handled by a computational version of the 3 formalism |2^] . 

In a general null evolution, data at the worldtube would be extracted from a Cauchy evolution . In this example, 
since we are not matching to a Cauchy evolution, we extract data at the worldtube from an analytical Cauchy solution. 
The characteristic evolution is carried out using a code described and calibrated in ||, transformed into an ingoing 
null code according to the procedure presented in || . 

Inside the grid there is a black hole whose interior is partially excised at an inner boundary, the "hole" , which is 
taken to be either a marginally trapped surface or the Q-boundary. We need to evolve only those grid points which 
are outside a discrete version of the hole. At the same time, we also need to allow the hole to move through the 
grid. In order to accomplish this we use a 3 dimensional mask function. Each grid point is assigned the value 1 by 
the mask if it either neighbors or is exterior to the boundary. All other grid points are masked to zero. The metric 
functions are evolved at each point with mask value 1 using data only from other points with a mask value 1 (i.e. 
points outside or neighboring the hole use data only from other points outside or neighboring the hole). In case they 
are needed (see below), values of the metric functions at grid points which are nearest neighbors just inside the hole 
are extrapolated radially inward using points exterior to them with mask 1. Other interior points are ignored. 

After all the metric functions have been evolved, we locate the hole at the new time. We then recompute the mask 
function and continue. If the boundary moves out, we simply throw away the data which we just evolved. If the 
boundary moves in, we have data at the new point which was obtained by extrapolation. Using this approach saves 
us from having to figure out if we have any points which were previously in the hole but are now outside. Such points 
automatically have data. It should be noted that we can safely assume that the black hole boundary will never move 
more than a grid point during any iteration. If it did, it would violate the Courant-Friedrichs-Lewy condition, which 
is built into the characteristic code to insure stability. 

The procedure for locating the trapping boundary is fairly simple. Basically, we use the previous position of the 
horizon (R) as a guess for its current location (if this is the first iteration, we use the position of the Q-boundary 
for the approximate location of the trapping boundary). Then, if we are finding the boundary using the parabolic 
relaxation technique based upon Eq.( |3.3| ), we compute 0; and then let R = R id — &iSX, and repeat until 1 1 z 1 1 2 is less 
than some threshold. If instead, we are using the minimization procedure, we compute F and then let R = R id — F5X, 
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and repeat until ||-F||2 is less than some threshold. The values of the stopping threshold and SX are parameters. The 
threshold can be safely set to be approximately Ar, (where Ar is the spacing between radial grid points) . The value 
of SX should be small enough so that FS\ < Ar but not so small that finding the horizon requires a large number of 
iterations. We have found it sufficient to choose a fixed SX at the start of the calculation, however it is possible to 
design a scheme in which SX is self-tuning and changes its value to speed convergence of each attempt to locate the 
horizon. 

A stability analysis of the explicit scheme (first order in time, second order in space) used to solve Eq. (3.3) shows 
that SX must scale as A£ 2 (where A£ is the spacing between angular points). This requirement would suggest that 
the proposed method is computationally expensive. However, our results show that the use of the horizon finder 
introduces a negligible overhead when dealing with long evolutions. With the aforementioned strategy of using the 
position of the Q-boundary as the initial guess, the finder may need many iterations to home in on the trapped 
surface, but it takes just a few iterations thereafter to track the surface. 

It is worth noting that, if instead of using the Q-boundary to determine the initial guess one uses a more educated 
guess for the location of the surface, the nu mber of steps can be dramatically decreased. This can be easily done in 



the case of the boosted black hole (see Sec. |V A| ). Using the expression for r obtained in that case, one can set as 



initial guess R — 2m/(cos0sinho! + cosh a), reducing the number of iterations on the first hypersurface from several 
hundred (using the Q-boundary) to less than 10 ( for va lues of a < 0.5). 



In the case of a wobbling black hole (see Sec. VD), we do not know the analytic expression for the marginally 



trapped surface even at the first hypersurface so that we use the Q-boundary as our initial guess. Table | shows the 
number of iterations made by the horizon finder for different values of the offset b and frequency u (for a characteristic 
grid having 45 x 21 2 points covering the space from r = to r — 4). After the first "find" the number of iterations 
necessary to track the hole is small; and since the finder solves a _/V| problem (as opposed to an N r x iV| problem), 
it does not add an appreciable computational time to the overall numerical evolution. 

In comparison, we obtain decid edly inferior efficiency in locating and tracking the hole by the minimization procedure 
using the flow given in Eq. ( 3.1 2| ) . Stability analysis of the finite difference method shows that SX must scale as A£ 4 in 
this case. Thus, although this minimization approach is attractive, it is not practical using finite difference techniques. 
It may be possible to improve the minimization algorithm by using pseudo-spectral techniques |n],^0| to calculate the 
flow but we have not explored this possibility in the present work. 



V. RESULTS 



Here we present some results of code runs for various initial conditions. We describe the physical behavior of the 
black hole in terms of the surface area of its marginally trapped surface. This surface area gives a measure of the 
energy of the radiation fields introduced in the initial null data. For the pure Schwarzschild case, the marginally 
trapped surfaces have area A s = IGttM^, in terms of the mass M s of the Schwarzschild black hole. More generally, 
the surface area of a marginally trapped surface determines its Hawking mass Mh |24| according to A = 16ttM^. 
Thus, on an ingoing null hypersurface Af Vl the function A(v) = M s — Mh(v) provides a measure of the energy between 
the marginally trapped surface and the worldtube. 

If a spacetime satisfies some suitable version of cosmic censorship, such as asymptotic predictability, and settles 
down to a Kerr black hole, then the area of any marginally trapped surface must be less than the area of the final 
black hole J2j|. In the present context, we do not have a global asymptotically flat solution so these results are not 
immediately applicable. However, if the black hole settles down to equilibrium with the Schwarzschild boundary 
condition on the worldtube, then at late advanced times we would expect A(v) — ► 0. 



A. Rotating Schwarzschild black holes 



Our first runs are for a Schwarzschild spaceti me described in coordinates rotating about the center of spherical 
symmetry, as described by the null data in Eq. (|]^). In this case, the evolved metric is known analytically and the 
marginally trapped surface is fixed at the horizon at r = 2m, so that convergence to the exact results can be checked. 
Our results confirm that the numerically evolved spacetime is accurate to second order in grid size. As expected, the 
horizon finder converges to the known location of the spherically symmetric marginally trapped surface. 
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B. Boosted Schwarzschild black holes 



A boosted Schwarzschild black hole provides a test of the ability to track the motion of a black hole and to calculate 
its surface area. Let x a be ingoing Eddington-Finkelstein coordinates. Define boosted coordinates x a by x = x, 
y — y, z — z cosh a — t sinh a and t = t cosh a — z sinh a. 

We locate an initial incoming null hypersurface Ao at t = — (x 2 + y 2 + z 2 ) 1 / 2 . Initial Schwarzschild null data 
corresponds to setting the Weyl data to zero at this initial time. Schwarzschild data at the extraction worldtube 
x 2 + y 2 + z 2 = R 2 is provided by transforming the metric to the x a coordinates. 

Let (v, f, 8, 4>) be standard ingoing null coordinates associated with the Cartesian Eddington-Finkelstein coordinates, 
and let (v, r, 8, <f>) be the null coordinates associated with the Cartesian coordinates x a by the extraction module. We 
synchronize them so that v — v — on Ao, which is the only null hypersurface common to the v and v foliations. 
Then the boost transformation implies that <fr = <fr and, on Ao, that 

cos 8 cosh a + sinh a 

cosv = — — — , (5.1) 

cosn a + cos 8 smii a 

with inverse 

cos 8 cosh a — sinh a 

cos0= — , (5.2) 

cosh a ~ cos 8 sinh a 

Calculation of the Jacobian of the angular transformation gives r = r (cos 8 sinh a + cosh a) — r / (cosh a — cos 8 sinh a) 
on the initial null hypersurface. 

In order to find the initial null data in the boosted null frame we must also relate v and v. Near Ao, we set v = 
(t + r) + 0(v 2 ) . Then by carrying out the transformation to leading order in v we obtain v = £>/(cosh a — sinh a cos 8) . 
This is enough to determine that initially J = and (3 = 0. 

To the next order, the null condition g a ^v iCc v t p = implies v = w/(cosha — sinh a cos 8) + kv 2 + 0(v 3 ) where 



sin 2 8 sinh 2 a 



(5.3) 



2f 2 (cosh a — sinh a cos 8) 3 

The extraction routine is based upon the gauge condition that v = t on the worldtube, so that v_t = 1 at r = R. This 



fixes the integration constant in Eq. (5.3) and gives 



sin 2 8 sinh 2 a 



2f (cosh a — sinh a cos 8) 3 
Similarly, the condition g ab v ya 8^ = that 8 be constant along the null rays implies 

cos 8 cosh a — sinh a 



(5.4) 



cosh a — cos 8 sinh a 



+ ~/v + 0(v 2 ), (5.5) 



where 



sinh a sin 8 

7,f = —■ (5-6) 

f 2 (cosh a — sinh a cos 8) 3 

Here the gauge condition buil t in to the extraction routine is that 8.t = on the worldtube. With this boundary 
condition, the integral of Eq. ( |5.6| ) gives 

sinh a sin 2 ^. 



f (cosh a — sinh a cos 8) 3 



The v dependence of r can now be obtained from the defining equation of a surface area coordinate r q = det(gAB), 
where q is the determinant of the unit sphere metric corresponding to the x A coordinates. This gives (r 4 q)^ = 
-r 4 qgABg AB ,<>, where 

g AB t , = 2x A , & x B M g &i . (5.8) 
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A straightforward calculation on the initial null cone v = gives 



■— — Q {ll,r + f 2 (cos0) §1§ ), (5.9) 
2 sin ' ' 



which, using Eq's. (5.5) and (5/7), reduces to 

r t $ = sinh a cos 6. (5.10) 

This determines the Jacobian of the transformation between the the stationary and boosted null frames at v = 0. 
Carrying out the transformation of the metric gives the initial null data for a boosted Schwarzschild black hole: 
J = j3 = U = and V = — r + 2ra(cosha + sinh a cos 0)~ 3 , where m is the Schwarzschild mass. 

After the hole has moved so that it is no longer centered about the vertex of the null cones, the null metric still 
has some simple properties at the poles (0 = and 9 = tt) due to the axisymmetry of the system; e.g. J = at the 
poles. This allows an analytic transformation between null and Kerr-Schild coordinates along the ingoing polar null 
geodesies ±z = —(t — T) + R, which lie on the null foliation and leave the worldtubc at t = T. Along these polar 
rays, v = T + R and the radial null coordinate is given by r = \z\ = — (i — T) + R. This allows us to reexpress the 
location of the poles of the horizon ±z = 2m analytically in null coordinates as 

2m ± v sinh a . 

T = 1 "T" • i 5.11) 

cosn a ± sinn a 

Since ti = on the initial null cone, the pole of the horizon hits the vertex of the null cone after retarded time 
v = 2m/ sinh a. 

As a test of the evolution and finder, the surface area of the boosted event horizon should remain constant. We 
observed that this is the case throughout the evolution, modulo the first order error introduced by the horizon finder. 
We have confirmed that the surface area converges to the value determined by the Schwarzschild mass as the grid 



spacing is refined. We have also checked that the poles of the horizon move in accord with (5.11), so that the pole 
which travels inward moves slightly faster than the pole moving outward. 

The algorithm performs the evolution and tracks the motion of the horizon stably, as long as the CFL condition 
is satisfied. Figure ^ shows a 2-D cut (at y = 0) of the horizon displaying the position of the hole at three different 
times. HJ 

C. Distorted black holes 

Here we study the approach to equilibrium of a distorted black hole with Schwarzschild boundary conditions on a 
stationary worldtube. First consider the case where the worldtube is static. We introduce a gravitational wave pulse, 
with compact support, on the first hypersurface by 



Ml—) (1--J \h^ L T2Y l , m ifre[R a ,R b ] 



J(v = 0,r,x A ) = { V T J \ r J V 2/ + 1 (5 12) 

otherwise, 

where 2^i,m is the spin two spherical harmonic, R a = 1.5, Rb = 3 and the amplitude factor A = 45. 

As the evolution proceeds, the pulse gets reflected by the outer boundary and eventually falls into the hole. Our 
results confirm the expected behavior of a black hole approaching equilibrium. Figure]^ shows the behavior of Mh(v). 
The surface area increases monotonically and approaches the value 16ttM 2 determined by the Schwarzschild mass of 
the exterior. 

We also introduced a pulse on the initial hypersurface in the case where the worldtube rotates (thus inducing a shift 
of its world lines with respect to the static Schwarzschild streamlines). We observed that at any given time, this does 
not result in any change in the location of the boundary. Hence, a rotating world tube does not affect the behavior 
of the Hawking mass confirming that the rotation is a pure gauge effect. 

D. A Wobbling Black Hole 

Beginning with the Schwarzschild metric in Kerr-Schild coordinates x a , we introduce the coordinates of an offset, 
rotating frame x a by t = t, x = (x + b) coswi — y sinwi, y = (x + b) s\nu)t + y cos tot and z = z. In this frame, we use 
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the metric and its derivatives on the world tube x 2 + y 2 + z 2 = R 2 to provide the boundary values for a characteristic 
initial value problem based upon a family of ingoing null hypersurfaces. Although the Schwarzschild metric is static, 
the worldtube wobbles relative to the static Killing trajectories. 

The initial value problem is completed by posing null data determined by setting the Weyl data to zero on the 
initial null hypersurface. On a non-spherically symmetric null hypersurface, the Schwarzschild Weyl data is no longer 
zero (and it is not possible to express it in simple analytical form) . Thus our choice of vanishing Weyl data introduces 
an amount of gravitational radiation which depends on the size of the offset. 

The resulting spacetime is neither spherically symmetric nor static. Relative to the worldtube, it describes a black 
hole wobbling and emitting gravitational radiation. Relative to the static Schwarzschild symmetry, the worldtube 
wobbles but the black hole still moves and radiates. This physical picture is confirmed in Sec. [v| by the behavior of 
the surface area of the marginally trapped surface. The results demonstrate that the region of spacetime interior to 
the world tube can be evolved stably by means of a characteristic evolution algorithm, when this interior region is 
truncated in the vicinity of a trapped region. This is illustrated in Fig. ^, which displays the maximum values of the 
norms of J and U over the entire grid vs. time. After an initial transient stage, they settle into a stationary state 
without any sign of instability whatsoever. 

Figure [| displays az = cut of the trapped surface at different times, showing the ability to track the movement 
of the hole by the horizon finder . As the evolution proceeds, the horizon "wobbles" through the computational grid 
with period T — 2tt/u>. We have evolved up to 2000M confirming this behavior. p(| 

The accuracy of the numerical evolution in the region exterior to the horizon is negligibly affected by the choice of 
using cither the Q-boundary or marginally trapped surface as the inner boundary. This is illustrated, for the wobbling 
case, in Fig. GL where we plot the values of J vs time at points outside the inner boundary, as obtained by both 
methods. The numerical values have a negligible difference. However, evolution with the Q-boundary is somewhat 
superior with respect to performance since no elliptic solver or other iteration procedure is required. 

The area of the marginally trapped surface again approaches equilibrium with the Schwarzschild exterior. This is 
illustrated in Fig. ^, where the surface monotnically increases and approaches a constant value (which converges to 
!6irM 2 in first order). The usefulness of the Hawking mass as a measure of energy is supported by the observation 
that A remains positive. 

VI. CONCLUSION 

Our results display many interesting aspects of black hole physics, although their physical understanding is not 
completely clear and would require a deeper study of the surface sources induced on the worldtube. The most 
important accomplishment of this work is that characteristic evolution is now ready to supply both the inner and outer 
boundary condition for the Cauchy evolution of black holes as soon as Cauchy-characteristic-matching is achieved. 
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APPENDIX A: SPIN- WEIGHTED EXPRESSIONS 

For the divergence of the outgoing rays, we obtain from equation ( |3.1[ ) that 

——9, = M{-V + r 2 (W + U >r BR) + r(e 2fi J/r)^{M) 2 - r(e 2/3 K /r), r (BR)BR - B[e 2f3 (KdR - JM)}}. (Al) 



For 'I', given in equation (3.11), we obtain 



* = ^{[-^(^ J ),- + ^e 2/3 J(e i e- 2/3 ), r ], r (§ J R) 2 + [^(QiK), r - ^ ] K{Qie~^) , r ], r (3ii)3i2 

- 2{e ie - 2l3 u ir ) tr dR - 2g(e, e - 2 ' 3 i/ r ) + [^e 2l3 {.m - xg)(e /e - 2/3 )], r §i? 

9 9 A p 2 @ 9 

+ 9[( — OiK), r BR - ( — O, J), r clR - {-<die- 2l3 ), r — {KM - JdR) + — e 2/3 (J3 - Kd)(& ie - 20 )] .} (A2) 
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TABLE I. Performance of the Horizon Finder 



lo b First hypersurface Second hypersurface After the 10th hypersurface 

1300 1 I 

0.05 0.05 1416 5 < 4 

0.1 0.1 1800 16 < 7 

0.2 0.2 1539 257 < 14 
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FIG. 1. J values at different locations obtained using the Q-boundary (Qb) or the Horizon finder (HF); the values obtained 
with the different methods proposed to excise the hole agree quite well, indicating not only that both methods can be applied 
for the excision but also that causality is respected. The characteristic grid had 45 radial points and 25 2 angular points while 
the Cauchy grid extended from —4 to 4 with 45 points in each direction. The offset was b = 0.2, the rotation frequency lo = 0.2 
and the mass of the Schwarzschild exterior was defined by m — 0.5. 
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FIG. 3. Tracking a wobbling black hole: Cuts of the horizon at z = at times (0.0, 11.0, 22.0) displaying how the horizon 
finder can track the movement of the black hole. The run was made with 45 radial grid points and 25 2 angular points. We also 
set m = 0.25, offset b = 0.1 and angular velocity w = 0.1. 
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FIG. 4. Stability of a wobbling black hole: Maximun values of \J\ and |?7| over the entire grid vs. time (for a run until 
v — 1400M) . After an initial stage, these values settle to a stationary state indicating the stability of the evolution. The run 
was made with 45 radial grid points and 25 2 angular points. We also set m = 0.25, offset b = 0.2 and angular velocity lo = 0.2. 
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FIG. 5. Behavior of the surface area vs. time for a distorted black hole: At late time the system approaches equilibrium. 
The amplitude of the pulse is A = 45 describing an I — 2, m = spin weight 2 pulse, extending from r = 1.5 to r — 3.0 at the 
first hypersurface. The code was run with 41 radial grid points and 17 2 angular points, setting m = 0.5 
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FIG. 6. Behavior of the surface area vs. time for a wobbling black hole: At late time, the area approaches a constant value. 
The run was made with 45 radial grid points and 21 2 angular points, with offset b — 0.1 and angular velocity lo = 0.1. The 
mass of the Schwarzschild exterior m was set to 0.5. 
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